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We derive a direct general map from the luminosity distance Dl(z) to the inhomogeneous matter distribution 
M(r) in the Lemaitre-Tolman-Bondi (LTB) cosmology and compute several examples. One of our examples 
explicitly demonstrates that it is possible to tune the LTB cosmological solution to approximately reproduce the 
luminosity distance curve of a flat FRW universe with a cosmological constant. We also discuss how smooth 
matter distributions can evolve into naked singularities due to shell crossing when the inhomogeneous "curva- 
ture" E(r) is a function which changes sign. 



I. INTRODUCTION 

Standard Friedmann-Robertson- Walker (FRW) cosmology is characterized by the following features: homogeneity 
and isotropy when averaged on 0(100) Mpc length scales, negligible spatial curvature, and a stress tensor which today 
has a chemical composition of approximately 70% dark energy and 30% pressureless dust. The standard inflationary 
embedding of FRW cosmology currently provides a successful picture of the universe. However, since the invocation 
of dark energy leads to a new coincidence problem, which is that dark energy dominance roughly coincides with the 
epoch of nonlinear structure formation, and since we still do not understand the cosmological constant problem, there 
has been a renewed interest in exploring alternative inhomogeneous cosmological models (which are not perturb atively 
related to FRW cosmologies) to see whether they might offer a competitively plausible picture of the universe GIHS 

iiisiiiaiii. 

Supernova luminosity distance measurements as a function of redshift offer compelling evidence for dark energy 
when interpreted in the context of FRW cosmologies t ll2lll3ll . To interpret such data in the context of inhomogeneous 
cosmologies, it is often not particularly useful to average the underlying inhomogeneous variables to obtain a forced 
interpretation in terms of FRW cosmology. This is because generically there is no preferred spatial slicing to compute 
averages and no meaningful map between the physical observables and the time derivatives of the spatially averaged 
variables. In some sense, the observables contain more information than that which could be carried by averaged 
variables, and hence determining the appropriate smearing map on the observables which matches the information 
that could be carried by averaged variables requires a knowledge of the underlying inhomogeneities in the absence of 
special symmetries. Therefore, to characterize inhomogeneous cosmologies, it is useful to compute observables such 
as the luminosity distance function directly in terms of variables describing the underlying inhomogeneous geometry. 

The Lemaitre-Tolman-Bondi (LTB) solution corresponds to a spherically symmetric exact solution to the Einstein 
equations with pressureless ideal fluid. Extensive analyses have been carried out for this model because it allows for 
investigations of inhomogeneities that cannot be represented as perturbative deviations from FRW cosmologies. The 
LTB solution is fixed by choosing three smooth functions {E(r), M(r), Ro(r)}, which allow an infinite number of 
different radial inhomogeneities. Most of the previous attempts to compute a cosmologically plausible luminosity 
distance Dr 1 (z) as a function of redshift z consisted of finding a map from {E(r), M(r), Ro(r)} to Dl(z) or a Taylor 
expansion of D L (z) about z = 00H10. 

In this paper, we derive a map from {E(r), Dl(z), Ro(r)} to M(r). This has the advantage that the observed 
luminosity distance function Dl(z) more directly dictates the underlying cosmological model, as opposed to having to 
guess the right {E(r), M (r), Ro(r)} to produce the desired Dl(z). In a similar inverse problem (with a different 
choice of variables) was considered which focused on the situation with E = 0, while in this paper, we will keep E 
general. In using this new map, we find the interesting fact that the luminosity distance function is typically (depending 
on the choice of E{r)) an effective probe of the LTB geometry only for z < 1, since the luminosity distance function 
has a universal behavior in the limit that {Rq = 0, E — > 0}. More precisely, Dr,(z) is fixed independently of M(r) 
for {i?o = 0, E = 0}, since that limit corresponds to the = 1 FRW universe. A negative consequence of this 
feature is that the differential equation map from {E(r), Dl(z), i?o(r)} to M(r) fails to be a numerically accurate 
map beyond z ~ 1. Nonetheless, we show how the numerical solution can be patched to a semi-analytic solution 
beyond z = 1 to obtain a good fit (to within around 5% for the redshifts of interest) to obtain an LTB cosmology 
which reproduces the luminosity distance function of an FRW cosmology with = 0.7 and ft^j = 0.3. 

Furthermore, we make a subsidiary observation about the class of LTB solutions which can mimic the observed 
Dl(z). We find that if the radial inhomogeneity profile goes from E{r) > to E(r) < as r increases while M'(r) 
is positive in that region, there is generically a danger of forming naked singularities, which can be interpreted as due 
to shell crossing. 
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The order of presentation will be as follows. In the next section, we review the conventional approach to obtaining 
the luminosity distance as a function of redshift in the LTB cosmologies. In Seclffl] we construct a set of differential 
equations (which we will refer to as the inversion method) which can be used to map the luminosity distance function 
into a particular LTB geometry. Afterwards, we apply the method to several examples. In Sec. [V] we discuss how 
smooth geometries can evolve into naked singularities when E(r) switches signs. We then summarize and conclude. 
For completeness, we present the FRW luminosity distance with Q,m + = 1 in Appendix A. In Appendix B, we 
write the LTB solution explicitly in a form that is not commonly found in the literature. 



II. CONVENTIONAL APPROACH 



The spherically symmetric Lemaitre-Tolman-Bondi (LTB) metric 



ds 2 = dt 2 - - ( R ' r l, , dr 2 ~ R 2 dn 2 (1) 
1 + 2E(r) 



satisfies the Einstein equation with 



Tt = Diagonal^ = -!±-AJ-~P = 0, -P = 0, -P = 0] (2) 



for differentiable functions E(r) and M(r) if R(t, r) satisfies the partial differential equation 

'd t R\ 2 _ 2E(r) 2M(r) 



R R 2 R 3 



(3) 



The function E(r) can be thought of as a generalized version of the spatial curvature parameter (E(r) oc ~kr 2 in 
FRW), while M[r) can be thought of as a generalized version of mass (for matter domination in FRW M(r) oc 
Pi<i 3 r 3 1 M 2 t , where pi is an initial energy density, a,; is an initial scale factor, and r is the radial coordinate). The 
function R in the FRW limit takes the form rait). The boundary condition for Eq. Q is provided by the radial 
function Rq (r) = R(to, r), where to is the time at which the boundary condition is set (note that this time is generically 
different from today). 

The luminosity distance in LTB model is approximately given l^ fHSll by 

D L (z) = (l + z) 2 R(t(z),r(z)) (4) 



dr = y/l + 2E(r(z)) 

dz (l + z)d t d r R(t(z),r(z)) { ' 

dt _ -\d r R(t(z),r(z))\ 



dz (l+z)d t d r R(t(z),r(z))' 



(6) 



where t(z) and r(z) physically represent the geodesic of the photon coming to us (located at r = 0) starting from the 
radial distance of our horizon. Any photon we observe that starts from a closer radial distance will have a redshift 
which is the same as that experienced by the horizon photon when the ratio of the frequencies measured from any 
closer radial position is accounted for since redshift is independent of the frequency. Given an arbitrary choice of 
M(r) and E(r), Eqs. (0 and (|6j are conventionally solved to obtain the luminosity distance function through Eq. 
Angular diameter distance can be obtained from the luminosity distance function by dividing by (1 + z) 2 . 

One can simplify Eqs. and in the regime in which <9 t P {i.e., the "local expansion" rate) maintains the same 
sign by rewriting Eq. Q as 



/ , . 2M(r) 

d * R = s f E ^ + W7y (7) 
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where s = ±1 specifies whether there is local expansion or contraction. The solution to this differential equation 
requires a specification of a function of r at the initial time hypersurface. We will define that function to be Ra(r): 
Ro(r) = R(to,r) (recall to is not necessarily today). Hence, we compute 

n 1 DM ^ f \\ E'{r)+M'/R-Md r R/R 2 

dtd r R(t(z),r(z)) = s = , (8) 

and rewrite Eqs. and as 



dr 



s y /l + 2E(r(z)) y j2E(r) + ^ 
dz ~ (1 + z)[E'{r) + M'/R - Md r R/R 2 ] 



(9) 



dt -s\d r R(t(z), r{z))\^2E{r) + -gL- 
dz ~ (1 + z)[E'{r) + M'/R - MdrR/R 2 } 



(10) 



These have the advantage that there are no second derivatives appearing in the equations, but have the added assump- 
tion that the sign s is a constant. 



IH. INVERSION METHOD 

In Sec.|n| we explained the conventional approach of obtaining the luminosity distance function Dl(z) for a given 
{M(r), E(r), Ro(r)}. In this section, we wish to stipulate Dl(z) and solve for the class of {M(r), E(r), i?o(r)} 
that corresponds to this luminosity distance. In particular, we will solve for M(r) for a given {E(r), Dl(z), Ro(r)}. 
This inversion method has the advantage that the physical observable Dl(z) can be mapped to the geometry of the 
underlying model directly without having to guess M(r). The physics is simply that if one knows a single radial 
geodesic history of a photon which was emitted at an event (ii,ri) and observed at (£2,^*2), one knows the full 
spacetime geometry in the region (t\ < t < f 2 , T\ < t < r 2 ) of the LTB solution owing to its spherical symmetry. 

The basic equations for this goal also stem from Eqs. (|9), dl0> . and the equation for luminosity distance, Eq. (0}: 

R(z) = R(t(z),r(z)) = -!^%. (11) 
(1 + zy 

Since Eqs. and dim depend on d r R, we would like to find an expression for d r R as a function of z. To this end, 
we use the exact solution, which is given in Eq. JB It : 



, , M(r) V -^(A r 

(t - t )V2E(r) + y/W(j~r)y/E{r)R(t, r) + M{r) - Q(r) = -±L ln[ 



M(r) 
E(r) 



■R(t,r) 



(12) 



Q(r) = y/Ro(r)y/E(r)Ro(r)+M{r) (13) 
(or Eq. JB3I if E < 0). Taking d r of Eq. d!2i . we obtain a linear equation for d r R, which can be solved to find 



o,n = t^R'jr) + ( 2V - 2{t - U) + SL-SZ-^ ^m+m 



f[Ro]VR 

M 2 

+ W 2 

M'(r) 

H — 

2E 



2Ey/R 
1 



f[Ro} f[R] E*/ 2 ^ER- + f[R } 
1 



M a + ER + y/ERHf [Rq} M 1 + ER + \f~ERj [R] 

r ER + f[R} 



(14) 



1 , nzm 



f[R }^R VER 



2 In 



Mi 



Mi 



VER^ + f[Ro] f 2 [R} + VERf[R} f 2 [Ro] + VER^f{Ro} 
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where 

f[X] = ^Mx + E(r(z))X (15) 

Mi{z) = M(r(z)). (16) 
Furthermore, the function M'(r) in Eqs. (J9jl and dlOl can be replaced by 

M'{r) = (17) 
dz dz 

Since there are three unknown functions {Mi(z), r(z),t(z)}, and Eqs. and (II Oi (with appropriate substitutions for 
9 r i? and M (r)) provide only two independent equations, we require another independent equation. This is provided 
by dR/dz through the chain rule: 



d / 2Mi dt „ dr 

— = si2E+—±— + d r R—. (18) 

az V K dz dz 

For a given set of {J5(r), Dl(z), Ro(r)}, the set of differential equations Eq. Eq. I ll Oi l, and Eq. dl 81 can be solved 
for {t(z),r(z), M\ (z)}. Finally, to obtain M (r), we invert r(z) to obtain 

M(r) =Mx{z{r)). (19) 

In practice, as we discuss below, the procedure we just described is a bit more difficult because the differential equation 
can become singular for certain choices of {E(r), Dl(z), Ro(r)}. Also, for numerical implementation, it is useful to 
write Eqs. (|9}, (fTOI . and ( fT8l in the form 



dl 



dz 



(A) (20) 



where (A) is a matrix that does not contain any derivative terms. However, this matrix contains hundreds of terms 
consisting of combinations of {E, E', Dl, ■4-Dl,M\,Rq}, and is not very illuminating in the general case. 

Regarding the initial conditions, note that because the differential equation also generically has a 0/0 division near 
r = 0, a numerical implementation must typically set the boundary condition at a small but nonvanishing r. For this 
purpose, it is useful to linearize the system about r = to obtain intuition about the boundary condition near the 
origin. If we assume that R(t, 0) = 0, E(r — 0) = 0, lim,.^ = 0, and R(t, r) w rd r R(t, 0) near r — 0, we find 

t-U = -z^\ z=0 + O(z 2 ) (21) 

where t{ is the value of t at z — 0. Unfortunately, similar expressions for M\ (z) and r(z) near z = depend upon the 
assumption of the scaling behavior of M(r) and E(r) near r = 0, and even the limiting expressions are algebraically 
complicated partly because of the presence of logs. For example, taking the ansatz E(r) oc r 2 , r = r'(0)z, and 
M\ (z) = M c z 3 , we obtain from Eq. (I20> a nonlinear consistency equation near z = for {r'(0), M c }. Hence, it is 
simpler to directly solve for the initial conditions numerically using Eq. (12 1 i . the ansatz r = r'(0)z, and Eq. ( 1201 in 
the limit z — > 0. 



IV. EXAMPLES 



In this section, we will solve Eqs. (|9}, (I10> . and J 1 81 for a variety of choices of {-B(r), Dl(z), Ro(r)} to demonstrate 
the inversion method described in Section||ll]in physically relevant examples. 
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A. Flat FRW example 

As a first example of mapping {Dl(z), E(r)} to {M(r), R(t,r)}, consider the matter dominated (A = 0) flat 
Friedmann-Robertson- Walker (FRW) universe 

{R Q (r) =0,E(r) = 0} (22) 
with s = 1. The function d r R can be found from Eq. dl4t as 

1 M'(r) 

W( ' ,r) = 3MW fl(t ' r) - (23) 
The resulting differential equations (Eqs. (0, JTOl . and PI ) are 



2 M 

dr V K(t,r) 



dz (l + z)[4r&M/R-±M'/R] 



y/2MiR dr 
(l + z)[§A Ml ]dI 

dt -R?l 2 



dz (l + z) % /2M7 



(24) 
(25) 

(26) 



dz 3 Mi V R dz 



Note that the dr/dz equation in this limit becomes independent of dr/dz ^ because of the absence of E. Hence, 
this equation can be written as 



Note that Eqs. (|27) and J28) give 



( 1+ ^^ hMi i=te' (28) 



M 1/3 

R = c^- (29) 



1 + 2 

,1/3 



where c=R{z = 0)/M{'°(z = 0). Solving for Mi and i?, we find 

Mi = 2 3/2 c 3/2 | 1 - 



1 + 2 



and 



"-ST^-VlTI 1 - (30) 

Remarkably, M(r(z)) = Mi(z) is fixed without specifying r(z). Note also that if we identify -\/2c 3 / 2 = -g-, this 
i?(z) corresponds to the function implied by the luminosity distance of an = 1 FRW universe. Indeed, inserting 
Eq. d23l > into the expression for Too in Eq- 0> one finds that the stress tensor corresponds to a homogeneous FRW 
universe: Too °t t~ 2 independently of r. Hence, this limit corresponds to the homogeneous matter density FRW limit. 



6 



If we choose 



r(z) oc (1 



1 



which corresponds to the geodesic -redshift relationship in an FRW universe, we obtain 

M(r) = Mi{z{r)) oc r 3 , 



(31) 



(32) 



which corresponds to the M(r) leading to the familiar FRW solution. From Eq. d!2i . we can solve for R finding 
R oc ar. However, note that we can choose another r(z) and obtain an infinite class of different M(r) functions 
corresponding to the same luminosity distance function implied by Eq. (1301 . 

We can turn this result around. As long as M/R ^> E, the luminosity distance curve no longer accurately probes 
the geometry of the LTB model since different geometries lead to approximately the same R(z) — Dl(z)/(1 + z) 2 . 
In particular, this means the inversion method will necessarily be unstable once the curvature term E can be neglected. 
Schematically, we will have 



dr 

dz 



Mi 



V2M77? dr 



(33) 



when ER/Mi becomes small and F ~ Since 



y/2MlR 



dr 

dz 



ER/Mi 



1 in the limit that ER/Mi — > 0, we have 



y/2Mi R 



F ~ - 





0' 



(34) 



which creates an unstable differential equation for 4| in this limit. 



B. Void Model 



We will now reproduce the void model of |4] by applying the inversion method as a nontrivial check. 
First, let us review the solution of |4]. Their ansatz for E{r) and M (r) can be written as 

^M-^i,o^ 2 (/3o-^[l-tanh^ii]) (35) 
M{r) = ^i )0 r 3 (a - ^[1 - tanh^]) (36) 



R(trec,r) = a rec r, (37) 

with 

{a = l,/3o = 0, A/3 = -Aq = -0.9, Ar = 0.4r , r « -=-, H x>0 « H }, (38) 

where Hq w 50km/s/Mpc and a rec ~ 10~ 3 is the effective scale factor at recombination. Note that in their solution, 
E(r) > never changes sign. Hence, the solution is always in a "locally open" universe (recall E(r) ~ — kr 2 in the 

FRW limit). The resulting luminosity distance Z)^ 01 °(^) fits the supernovae luminosity distance well 0]. 

To use the inversion method to derive this M(r), we set E(r) to Eq. ( I35i . -R(z) = D^°^ c '(z)/(1 + z) 2 , and 
Ro(r) equal to Eq. ( 1371 and solve the differential equations Eqs. (|9}, il 0> . and dl 81 subject to the boundary condition 
« 0.855^-, r(zj) w 0, Mi(z,) w 0}, where ^ sa and lim r ^ Mi(z(r))/r 3 sa 0.084^ was taken to match 
Eq. J38i [ 19]. These boundary conditions correspond to ray tracing starting from the "center of the universe" where 
observations are assumed to be made. 

The results of the r(z) and M(r) reconstruction are shown in Fig.^ The solid curve was constructed using the 
inversion method and should ideally match the dashed curve. However, the inversion method breaks down for z ~ 0.6 
because there ER/M — > 0, in which case the dr/dz equation becomes unstable as explained in Subsection llV Al To 
see this explicitly, we plot ER/M as a function of redshift z in Fig.|2] 
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Figure 1: In the graph on the left, the solid curve shows the reconstructed r(z) (in units of 1/Ho) through the inversion method and 
the dashed curve shows the r(z) that would be obtained if M(r) of Eq. I36i were given as an input. The breakdown of the inversion 
method for this model around z ss 0.6 is apparent and is explained in the text. In the graph on the right, the solid curve gives the 
M (r) as determined through the numerical evaluation of Mi(z(r)), and the dashed curve gives M (r) as given by Eq. 1361 . 




Figure 2: ER/M versus photon redshift z. For z = 0.6, ER/M ~ 10 2 and the dr/dz differential equation is unstable. Another 
way to view the instability is that the luminosity distance function is not sensitive to the geometry for z > 0.6 in this model. 



C. Cosmological Constant without Cosmological Constant 

Until now, we have given two examples of using the inversion method to reproduce known LTB models. In this 
subsection, we construct an LTB model which reproduces an identical luminosity distance to that produced by an FRW 
universe with Qm = 0.3 and ft a = 0.7 up to a finite redshift. In particular, we will set R(z) = Dl(z)/(1 + z) 2 , with 
Dl{z) given by Eq. JAli . 

Although the inversion method is unstable for z > z c , where z c is a cutoff redshift (whose generic existence is 
suggested by the two examples previously presented), we can use the inversion method for z < z c and then smoothly 
patch that model (i.e. the function M(r)) on to a function M{r) ~ r 3 . Since the luminosity distance function is 
sensitive to the cosmological constant only for z < 1, such models tend to give good fits to the observed luminosity 
distance data. 

For our trial example, in addition to setting R(z) according to Eq. \A\\ and Ro(r) = 0, we choose 

E(r) = itf Vexp[-2tf r], 

which has a local spatial curvature only for r rts 0. The solution t{z), r(z), M%(z) is obtained by solving Eqs. (|9j, 
(I10> . and (1181 with the initial conditions to the differential equations set as outlined near Eq. \2\\ . Explicitly, we take 
t{zi) = 0.93/H , r(zi) = 0.97zi/H , M t (zi) = 0.028r 3 (^), and z, = 10~ 5 . As previously discussed near Eq. i2lt . 
we solved linearized equations to obtain t{zi) and r(zi), assuming M\{zi) = 0.028r 3 (zi). (The coefficient 0.28 in 
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Figure 3: Using the inversion method, we have computed the energy density p(t, r) corresponding to the luminosity distance of an 
FRW universe with Qm ~ 0.3 and Qa = 0.7. The solid, dashed, and dotted curves correspond to p(t, r) at times t to( j a y, ^today/^' 
and ^today/3' respectively. With Ho — 70 km/s/Mpc, ijoday * s ^ billion years. The redshift corresponding to r ~ 0.5/Ho is 
z ~ 0.4. The energy density at r = is nonzero: 47rp(t t0[ j a y, 0) / (M^Hq) ~ 10 -3 . 

. 4 




-0 . l l 

1 2 3 4 5 

z 



Figure 4: The solid curve corresponds to (D^ 0,7 (z) — D m (z)) / D^ ' 7 (z) which corresponds to the deviation of the numer- 
ically constructed "model" luminosity distance curve from the FRW curve (where the FRW model is one with {Qa = 0.7, Qm — 
0.3}). For comparison, we plot the dashed curve, which corresponds to (D^°' r (z) — (z)) / D^ ' 7 (z) , where D£~ cor- 

responds to a flat matter-dominated FRW model without any cosmological constant. It is clear that the model reproduces the 
luminosity curve exactly from z — up to z = 0.4 (by construction) and there is a less than about 5% deviation from the FRW 
{Qa = 0.7, Qm = 0.3} luminosity distance curve until about z = 3. Furthermore, the deviation error is seen to plateau at large 
redshifts. 



Mi(zi) is related to the energy density at r = through Eq. (0.) The differential equation was then solved from 
z w to z — z c w 0.4. From approximately this redshift onward 1 20], the inversion method is unstable because of the 
form 0/0. (It is unclear from our analysis whether this instability is purely numerical or whether this instability can be 
removed by a reformulation of the equation and perturbations of the boundary conditions. We defer this question to a 
future work.) Hence, we fit a smooth M(r) function starting from close to this point: 

_ f M^ m (z num (r)) r < 0.5/H 
M(r) ~ { ^(r + Cs) 3 r>0.5/H 

where the superscript "num" refers to the numerical solutions to Eqs. l|9}, i ll Oi l, and Jl 8i and C\ and C2 are adjusted to 
match the value and the r derivative at r — 0.5. We then take the extended M(r) and solve Eq. (IB 1 1 to obtain R(t, r) 
numerically. Finally, we then solve Eqs. @ and l|5} to obtain the full luminosity distance beyond z = z c w 0.4. Note 
that a functional choice of C\ [r + C2) 3 was made to obtain an approximately homogeneous p(t, r) for r ~ 1 /Hq. 

The resulting solutions can be seen in Figs.[3]and|4] Clearly, the inhomogeneous model, whose energy density is 
shown in Fig0 reproduces the luminosity distance of the {Qa — 0.7, f2 m — 0.3} FRW model exactly from z = to 



9 



z = 0.4 and the luminosity distance curve deviation from that of the FRW model is less than around 5% until z = 3. 
(The {Ha = 0.7, flj^ = 0.3} FRW model is known to give a good fit to the supernova data.) 

Hence, we have explicitly demonstrated that the LTB model can be tuned to obtain a luminosity-distance-redshift 
relationship which accurately reproduces that of a standard flat FRW universe with a cosmological constant and dark 
matter. This solution differs from previously proposed solutions in that the luminosity distance curve is exactly that of 
the FRW universe with a cosmological constant from z = to z — 0.4. It is interesting to note that the energy density 
is approximately homogeneous for large r but has a void close to r = (see Fig. O, similarly to the model of |4|. 



V. NAKED SINGULARITY FORMATION 



In this section, we demonstrate that the LTB model is susceptible to the formation of naked singularities when E(r) 
switches sign from positive to negative as r increases while M'(r) is positive in that region. More specifically, we 
consider situations in which k(r) ~ —2E(r)/r 2 (i.e., the "local" spatial curvature factor) makes a smooth transition 
from an "underdense" universe to an "overdense" universe as r crosses r from below. This naked singularity may 
be interpreted as due to formation of caustics arising from matter accretion. Since any realistic system has nonzero 
pressure at sufficiently small length scales, this singularity is unphysical and is an artifact of pressureless dust approx- 
imation. This type of naked singularity would develop if one naively smooths out inhomogeneity profiles such as the 
one used by |2|]. 

Now, note that because of the expression for the energy density in Eq. (|2j, for the energy density to be positive, 
M'(r) must have the same sign as R, r . Using Eq. iB5\ . which is valid for ER/M <C 1 (which is almost always 
true in the vicinity of r = ro, the point at which E(r) changes sign smoothly), let us consider the situation in which 
Ro( r ) = and t = (the case considered by [2]). We find 



R, r = P(t,r) 



E'(r) - 



E(r)M'{r) 5 • 2 1 / 3 M'(r) 



3 M(r) 3 ■ 3 2 /3 £2/3 M i/3( r ) 



(39) 



where 



P(t,r) = 



3 ■ 3 1 / 3 f 4 / 3 
5 . 2 2/3 MV3( r ) 



(40) 



Supposing that the last term proportional to M' (r) / (i 2 / 3 M 1 / 3 ) can be neglected compared to the E'{r) term and 
M'(r) > 0, we have the condition 



E'{r) > 



E(r) M'(r) 
3 M(r) 



(41) 



if M > 0. If we choose a smooth E(r) which goes from a positive value to a negative value as r increases past ro, the 
Taylor expansion of E(r) is then required to have the form 



E{r) = Bi(r )(r - r ) + ^E 2 (r )(r - r f + 



(42) 



with Ei < 0, where E n corresponds to the nth derivative of E(r). Hence, in a sufficiently small neighborhood of 
r = r , for Ei (r ) ^ we have the condition 



(r-r )M'(r) 
3 M(r) 



> 1. 



(43) 



If we assume M'/M does not switch sign at r = ro, this condition is clearly not satisfied in the vicinity of ro. This 
means that when the E'(r) term dominates the R, r expression, the energy density cannot remain positive definite. 

Now, suppose that the M'(r) / (t 2 ^ 3 M 1//3 ) term dominates over the E' term in Eq. (1391 . which would certainly be 
true near the big43ang singularity at t — > 0. We find the following condition for M' > and M > 



5 ■ 2 1 / 3 M 2 / 3 (r) 
3 2 / 3 t 2 l 3 



> E{r) 



(44) 
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Figure 5: The energy density is plotted as a function of radius for times t = 0.01L (solid) and t = 0.042L (dotted) for the model 
of Eqs. 1451 and 1461 with {ro = 0.7L, k = 50, po = 3.5-^}. As expected, the energy density diverges as R, r — > at r ~ ro 
near the time t > 0.04L. Note also that the energy density for r > ro is larger than r < ro, since one side is overdense while the 
other side is underdense. 



which can be satisfied in the vicinity of r = tq. This means that energy can be positive definite in these regime. 
Hence, we arrive at a naively puzzling question why the energy density, which is initially positive definite near the big 
bang singularity, develops into a negative energy density when E'(r) term governs the value of p. The answer is that 
a naked singularity develops during the course of matter evolution. 
To see this in a more obvious way, let us take a concrete model of 

~3~3 

M(r)=/*^- (45) 
6 

and 

E{r) =2^ tanh ( K ^-^ ) . (46) 

which represents a version of the model of |2] with the step function smoothed out. We can explicitly solve for the 
coordinate where R, r becomes negative by solving R, r = : 

3 .' 6 V, L 'fT ~ - «anh[^pl] - 0. (47, 

Expanding about r = r to quadratic order in r — r , we find the radius at which R, r becomes zero to be 




'" ' ' / - 5 " 6a/BX -[ a g^ /8 *-»/. + -^ ff + l. (48) 



Requiring that the solution be real results in the condition 

*>V5[poa^ 3 ]( — ) 3/2 , (49) 
Kr 

T 2 

after which time, R, r has a value of zero near r ~ ro H — . Comparing this time with the time t c at which the 



K^ro 
6~ 



overdense part of the universe starts to collapse (£ c sw """^ Po ), we see that if k ^> 1 while O(ro) ~ O(L), i?, r will 
reach zero before t reaches t c . 

When R, r reaches zero at nonzero r = r s without M'(r) vanishing, we have a Ricci curvature singularity there. We 
have also checked the divergence of the energy density near the spacetime region of the naked singularity for several 
numerical examples, one of which is shown in Fig. [5] Since no metric element became zero before the singularity 
comes into existence, we see that the curvature singularity is not protected by a horizon; i.e., it is a naked singularity. 
As is well known 117T . naked singularities can develop with perfect fluid systems. Hence, if E(r) changes sign at 
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Figure 6: The energy density as a function of radius for three different times. The dotted curve corresponds to t = (recombina- 
tion), the dashed curve corresponds to t = 0.1/Hq, and the solid curve corresponds to t — 0.5/Hq. 



r = r with M{r) > and M'(r) > in the vicinity of r , a naked singularity develops, and the Einstein equations 
break down in that region. Furthermore, if we naively continue using the Einstein equations past the development of 
naked singularities, one finds a negative energy density region which is unphysical as expected. 

Here we note that the model of |4] does not have any naked singularity problems associated with shell crossing. To 
see this, we will use the analytic approximation of Eq. JB5I applicable to the region in which \E(r)R/M\ <C 1. In the 
scenario specified by Eqs. i38l , we have 

\ E(r °M^ r0) \ « |0-3(10- 4 + AH tf'\l - J*™^ + 0.06[10- 4 + 47^/3)1 « i (50 ) 

until the time t > 1/Hq. One can then use Eq. JB5I in this regime to evaluate R, r jr and find that the energy density 
does not diverge anywhere. The key reason is that E(r) never switches sign, unlike the previous case. 

There is a peculiar physical feature of the energy density as a function of time, as shown in Figure |6] Some of the 
matter density from the transition region near tq initially transfers to the r = region before being diluted away by 
curvature dominated expansion. In other words, the density near r = is not a monotonic function of time because 
the initial conditions are set up such that the fluid is flowing towards r = 0. 



VI. CONCLUSION 



In this paper, we have derived a set of differential equations in the context of LTB cosmologies (spherically sym- 
metric pressureless dust solutions to the Einstein equations) which can be numerically solved to obtain almost any 
luminosity-distance-redshift relationship that can be produced by a homogeneous and isotropic FRW model. We have 
solved this set of differential equations for several examples to demonstrate the feasibility of our inversion method. 
Unlike many other methods in the literature, our method can be used to dial in the geometry that generates the desired 
luminosity distance exactly in a finite redshift interval. We have also given explicit examples of naked singularity 
formation in LTB cosmologies in the region where E(r) changes sign. 

Although our work serves as a step towards building cosmological models competitive to the standard FRW cos- 
mology, the research program is far from completion. Toy models such as the LTB cosmologies are arguably not 
yet convincing contenders for compelling alternatives to standard inflationary cosmology, since being at the center 
of the universe requires giving up the Copernican principle (though recently there has been some effort to alleviate 
this problem ||7|]) and there is not yet a convincing structure formation history that could explain such radial inhomo- 
geneities. Nonetheless, it is not obvious whether this class of LTB models can be ruled out from current observations 
0. Furthermore, there is some evidence that our galaxy is in a void, qualitatively similar to the voids presented in two 
of our numerical examples. However, the existence of the void is still currently being investigated (see for example 
[18]), and whether the magnitude of the void is sufficient to cause the redshift effects presented in this paper is unclear. 
Regarding tests of LTB cosmologies, any physical probe testing the isotropy of the universe from a point separated 
sufficiently far away from us could in principle be useful. We leave this problem for future investigation. 
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Appendix A: FRW LUMINOSITY DISTANCE 

The FRW luminosity distance in flat FRW universe with pressureless dust fraction £1 m and cosmological constant 
fraction £!a without spatial curvature (11m + = 1) is given by 



l + z 



D L {z) = -jr- / dyl(y) (Al) 

-no Ji 

i(y) = 1 (A2) 

^n A + (l - n A )y* 

where the integral is numerically a value of order 1 with a logarithmic dependence on z. Note that a Taylor expansion 
is not very efficient approximation of the integral as fifth order expansion gives a good fit to only about z « 1.5. 
To check the plausibility of this expression, consider the solvable case of matter dominated universe Qj\ = 0: 

D L (z) = ±±£ [^dyy- 3 ^ (A3) 

"0 Ji 

= 2^-[l--=L=] (A4) 

which is the familiar expression. If this is matched to LTB model, we would write 

Dl(z) 
[l + zf 
1 , 1 1 

(A6) 



'H l l + z (l + z) 3 / 21 
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Appendix B: EXPLICIT SOLUTION 

Einstein equation Q can be solved by constructing characteristic curves. The result can be expressed as 



M (r) , yW^) + ^J^ + R(t,r) 



(t t )V2E(r) + ^Wy)VE(r)R(t,r)+M(r) - Q(r) = -j^L ln[ V r ] (Bl) 

Q(r) ee y/R (r)y/E(r)R (r)+M(r) (B2) 

for E > while for E 1 < 0, we have 



F(r) - (t - t )V2E(r) = ^= arcsinU -f^ y/R(t,r)] - y/R(t,r)y/E(r)R(t,r) + M(r) (B3) 



F(r) = arcsin[^ 1 ^ v^W] - ^i^/W^oM + M(r). (B4) 

Here, Ro(r) is the initial condition specification for R(t,r): i.e. R(to,r) = Ro(r). Note that the solution can be 
written a little more explicitly in the limit \EM/R\ < 1 and |_BM/ii| > 1. For |i?M/i?| < 1, we find 



. 2.2 1 /3i?3 + 6 .2 5 / 6 VMi?o /2 (<-to) + 9-2 1 /3 M (t-io) 2 2^ /2 

" = ; — ; tr, ; — 

5[2Rq /2 + 3V2M(t - t )} 4 / 3 5[2Rq' 2 + 3V2M{t - i )] 

For ER/M ^> 1, the leading order self-consistency requires 

R = R + V2E(t-to) + ^ln[^-} + 0(^). (B7) 
Hence, we can approximate in the regime of interest 



R^^2^)(t-t Q )+R (r) + ^\nll + ^^^(t-t )\ . (B8) 



